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ABSTRACT 

Giant molecular clouds (GMCs), where most stars form, may originate 
from self-gravitating instabilities in the interstellar medium. Using local three- 
dimensional magnetohydrodynamic simulations, we investigate ways in which 
galactic turbulence associated with the magnetorotational instability (MRI) may 
influence the formation and properties of these massive, self-gravitating clouds. 
Our disk models are vertically stratified with both gaseous and stellar grav- 
ity, and subject to uniform shear corresponding to a flat rotation curve. Initial 
magnetic flelds are assumed to be weak and purely vertical. For simplicity, we 
adopt an isothermal equation of state with sound speed Cg = 7kms~^. We flnd 
that MRI-driven turbulence develops rapidly, with the saturated-state Shakura 
& Sunyaev parameter a ~ (0.15 — 0.3) dominated by Maxwell stresses. Many of 
the dimensionless characteristics of the turbulence (e.g. the ratio of the Maxwell 
to Reynolds stresses) are similar to results from previous MRI studies of accre- 
tion disks, hence insensitive to the degree of vertical disk compression, shear 
rate, and the presence of self-gravity - although self-gravity enhances fluctuation 
amplitudes slightly. The density-weighted velocity dispersions in non- or weakly 
self-gravitating disks are ax ^ <7y (0.4 — 0.6)cs and (Xz ~ (0.2 — 0.3)cs, suggesting 
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that MRI can contribute significantly to the observed level of galactic turbulence. 
The saturated-state magnetic field strength B ~ 2/xG is similar to typical galactic 
values. When self-gravity is strong enough, MRI-driven high-amplitude density 
perturbations arc swing-amplified to form Jeans-mass (~ 10^ Mq) bound clouds. 
Compared to previous unmagnetized or strongly-magnetized disk models, the 
threshold for nonlinear instability in the present models occurs for surface den- 
sities at least 50% lower, corresponding to the Toomre parameter Qth ~ 1-6. We 
present evidence that self-gravitating clouds hke GMCs formed under conditions 
similar to our models can lose much of their original spin angular momenta by 
magnetic braking, preferentially via fields threading near-perpendicularly to their 
spin axes. Finally, wc discuss the present results within the larger theoretical and 
observational context, outlining directions for future study. 

Subject headings: galaxies: ISM — instabilities — ISM: kinematics and dynamics 
— ISM: magnetic fields — MHD — stars: formation 

1. Introduction 

It has long been a challenge to understand the origins of giant molecular clouds (GMCs). 
Observations and theory both favor "top-down" (instability) rather than "bottom-up" (co- 
agulation) mechanisms, at least for the most massive clouds (~ 10^ - 10*^ Mq) in which the 
majority of the molecular material is found and the preponderance of stars are born (e.g., 
Bhtz 1993). The principal difficulty with stochastic coagulation is that it proceeds much 
too slowly to achieve GMC masses given other constraints (e.g., Bhtz & Shu 1980); recent 
support for this view includes the observed lack of raw material in small clouds within spiral 
arms (Heyer & Terebey 1998), and the conclusion from numerical simulations that magne- 
tohydro dynamic (MHD) turbulence cannot prevent gravitational collapse within clouds that 
exceed the magnetic critical mass^ (e.g., Ostriker, Stone, & Gammic 2001; Heitsch, Mac 
Low, & Klessen 2001). The latter result implies that intermediate products of coagulation 
^ 10'^ Mq, if they were created, would undergo rapid star formation and be destroyed by 
the associated energy feedback before subsequent mergers could attain ~ 10^ M©. 

The two main classes of instabilities proposed for GMC formation are Parker and self- 
gravitating modes, and their hybrids (e.g., Elmegrcen 1995). Although the azimuthal wave- 
lengths and initial growth rates of Parker modes suggest they could prompt cloud formation 
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in spiral arms (e.g., Mouschovias, Shu, & Woodward 1974), numerical simulations have re- 
cently shown that nonlinear development favors short radial wavelengths and saturates at 
moderate amplitudes. The low-density-contrast structures that develop do not appear to 
be the primary progenitors of gravitationally-bound GMCs (J.Kim et al. 2000, 2001; Kim, 
Ostriker, & Stone 2002, hereafter Paper III). 

Direct sclf-gravitating instabilities have similar preferred azimuthal scales and growth 
rates to Parker modes, but arc not self-limiting in the nonlinear regime. In regions without 
spiral structure, epicyclic motion can temporarily conspire with shear in the swing amplifier 
process (Toomrc 1981) to yield gravitational runaway provided a sufficient surface density. 
Based on two-dimensional simulations, Kim & Ostriker (2001, hereafter Paper I) found a 
nonlinear instabihty threshold at Toomre Qth = 1-2 — 1.4 (for definition, see eq. [4] below) 
for a range of magnetic field strengths; from three-dimensional simulations with zero or 
thermal-equipartition magnetic fields. Paper III found Qth ~ 1- 

In high-density portions of spiral arms, local reduction or reversal of shear suppresses 
swing amphfication, but magneto- Jeans instability (see Lynden-Bell 1966; Elmegreen 1987 
and Paper I) can yield self-gravitating growth since magnetic torques between azimuthally- 
adjacent regions limit the epicyclic motion and spin-up of contracting condensations. Numer- 
ical simulations in thin-disk models (Kim & Ostriker 2002, hereafter Paper II) have shown 
that the magneto- Jeans instability in spiral arms first leads to growth of gaseous spurs, which 
then fragment into Jeans-scale condensations. 

The investigations of cloud formation in Papers I-III considered either two-dimensional 
or three-dimensional models with zero or strong magnetic fields (to study Parker modes). 
Under more realistic galactic conditions with subthermal mean magnetic fields, however, new 
dynamical effects are introduced by the action of the magnetorotational instability (MRl). 
Although the MRl has primarily been studied in the context of accretion disks, Sellwood and 
Balbus (1999) proposed that MRl and related processes may be important for maintaining 
turbulence in the outer H I disks of galaxies where the star formation rate (and supernova 
rate) is low. The recent calculations of Wolfire et al. (2003) showing that a cold phase 
must be present in the far outer Milky Way support the notion that MRI-driven turbulence 
prevents the settling and gravitational fragmentation of a thin, cold gas layer that would 
otherwise ensue. 

To understand why the MRl is likely to be important in galaxies, it is helpful to consider 
the instantaneous instability criterion for MRl modes oc e*^'^, 

(k-VA)^|k|2 dlnn 
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(e.g., Balbus & Hawley 1998). Here, 0,{R) is the disk's angular rotation profile, and 
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Va = B/(47rp)^/^ in terms of the local density p and mean magnetic field B. Physically, 
this criterion represents the need for (moderate) magnetic tension torques between radially- 
displaced fluid elements that are azimuthally sheared by the background flow to reinforce the 
initial perturbations in angular momentum (Balbus & Hawley 1991). In a galaxy without 
strong spiral structure or in interarm regions when a prominent spiral pattern exists, q ^ 1 
and unstable MRl modes with a variety of wavcnumbers will be present.^ 

The detailed dynamical characteristics and consequences of MRI-driven turbulence in 
galaxies have not previously been studied directly. Although the literature contains many 
theoretical and numerical treatments of MRI (e.g. Balbus & Hawley 1991; Hawley, Gammie, 
& Balbus 1995; Stone et al. 1996; Miller & Stone 2000; Hawley & Krolik 2002, see also Balbus 
& Hawley 1998; Stone et al. 2000 and references therein), these works have focused primarily 
on angular momentum transport in accretion disks around stars and black holes. These 
accretion disks differ significantly from galactic disks in that the latter are self-gravitating, 
have flat rotation curves rather than Keplerian rotation, and are more vertically compressed 
by both gaseous self-gravity and external stellar gravity rather than solely by a tidal stellar 
potential. How do the physical properties of a disk affect the level and characteristics of the 
saturated-state turbulence driven by MRl? How does this turbulence nonlinearly interact 
with self-gravity to form GMC-like structures in galactic disks? In this paper, we address 
these and related questions by performing and analyzing a suite of three-dimensional MHD 
simulations. 

As in Paper 111, our disk models arc local, isothermal, and vertically stratified. Initial 
magnetic fields are assumed to have a weak, uniform strength, pointing in the vertical direc- 
tion (but note that after MRI develops, B becomes toroidally-dominated) . The background 
potential corresponds to that yielding a flat galactic rotation curve. In this first work on 
galactic MRI, we do not consider other features such as stellar spiral structure and effects of 
the multi-phase ISM that would be present in real galaxies; these effects will be considered 
in subsequent work. 

In §2, we briefly describe our model parameters and computational methods. In §3.1, we 
investigate the characteristics of MRI-driven turbulence in our non- or weakly self-gravitating 
models, and compare with results from previous accretion disk simulations. We directly 
measure the density-weighted velocity dispersions to assess whether MRI can be a significant 



^E.G. using Vl = 26knis ^kpc ^ and scale height H ^ 150 pc, with primarily-toroidal mean Bq ^ 1.4/iG, 



and assuming wavelengths < A;?, Ag, eq. (1) requires 
azimuthal wavelengths » H and vertical wavelengths up to 
too strong). 



0.1 - which essentially calls for 
H (provided the vertical magnetic field is not 
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source of turbulence in the diffuse ISM. We also discuss the role played by self-gravity in 
supplying turbulent energy to the ISM. In §3.2, wc present results on Q thresholds for the 
formation of gravitationally bound clouds in MRI-unstable disks, with comparison to our 
model results from Paper III. The properties of the bound clouds formed and the role of 
magnetic fields in braking their rotational spins is briefly discussed as well. Finally, in 
§4 we conclude and discuss our findings and related issues within a larger context, noting 
outstanding questions related to galactic turbulence and GMC formation that will be the 
focus of future work. 



2. Methods and Model Pcirameters 

For local, three-dimensional simulations of vertically stratified, MRI-unstable, galactic 
gas disks, we use the same numerical code and model formulation described in Paper III, with 
the exception that we modify the strength and spatial distribution of the initial magnetic 
field. We adopt an isothermal condition in both space and time, with an effective isothermal 
speed of sound Cg — 1 kms~^, and solve for an initial vertical density distribution Pq{z) in 
hydrostatic equilibrium between the thermal pressure gradient and the total (gaseous plus 
stellar) gravity. For the strength of external stellar gravity, we define 

where c^^^ and are the vertical velocity dispersion and surface density of stars, respec- 
tively, while So denotes the gas surface density (see Paper III). To represent average disk 
conditions, we take sq = 1 corresponding to a^^z ~ 20 kms~\ ^ 35 M© pc~^, and 
Eo 13 Mq pc-2 (e.g., Kuijken & Gilmore 1989; Holmberg & Flynn 2000). Note that these 
values give comparable stellar-disk and gas-disk vertical gravity at \z\ ~ if, where H stands 
for the gas disk scale height (Paper III). 

We introduce an initial uniform magnetic field, Bq = Sqz, that points in the vertical 
(z) direction. We parameterize its strength using 

. I = 1^^^. (3) 

where vx = Bq/ [AT^p^y/'^ is the Alfven speed. In the numerical models presented in this 

paper, we adopt /9(0) = 100 or 400. The physical value of the initial midplanc magnetic field 
is thus ~ 0.2 — 0.4/iG. Since the density drops off with height, high-latitude regions may have 
(3{z) < 1 with which the growth of MRI is ineffective, but the weak midplane magnetic fields 
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guarantee that the most unstable wavelength of MRI is confined within \z\ ^ H. Note that 
although the initial magnetic field in our models is primarily vertical (a numerical expedient 
that yields rapid growth of MRI), this condition is not inconsistent with observations of pref- 
erentially azimuthal galactic B-fields, because the saturated-state MRI has mainly toroidal 
B, as we shall show in §3.1. Also note that the amplitudes of saturated-state turbulence 
depend on the net fiux for pure- 2; initial 5- fields (Hawley, Gammie, &; Balbus 1995, hereafter 
HGB); we will briefly discuss this in §3.1.1 and §4.1. Other classes of initial magnetic held 
distributions that have been studied in non-self-gravitating MRI simulations (e.g., HGB; 
Stone et al. 1996; Miller & Stone 2000) contain zero net magnetic flux. Since simulations 
which begin with zero net flux take much longer to reach saturated turbulence, we have fo- 
cused on uniform vertical field models in this paper. Further exploration of intermediate field 
strength and/or orientation may yield interesting results; we shall, however, defer broader 
parameter surveys to future investigations. 

An important measure of susceptibility to gravitational instability is the Toomre Q 
parameter (Toomre 1964) 

^ ^ ^ 4 / C_s \ / « \ ( So 

^ ttGEo ■ V7.0kms-V V36kms-ikpc-V Vl3 M0PC-2 

where k, is the epicyclic frequency. We evolve models with Q in the range of 1 to 2 to 
represent various conditions in disk galaxies (e.g., Martin & Kennicutt 2001). Note that the 
often-cited critical value of Qth = 1 applies to axisymmetric instability in an infinitesimally 
thin disk. Even in unmagnetized cases, finite disk thickness generally reduces the critical 
Q value for axisymmetric instability to Qth = 0.676 for purely self-gravitating disk with 
So = 00 (Goldreich & Lynden-Bell 1965a; Gammie 2001) and Qth — 0.75 for a disk with 
So = 1 (Paper III). Therefore, our model disks having Q ^ 1 and Sq = 1 are all stable to 
axisymmetric perturbations. 

We integrate the time-dependent, compressible, ideal MHD equations in a local Carte- 
sian reference frame whose center orbits the galaxy with a fixed angular velocity Q = ^{Ro) 
at a galactocentric radius Rq. In this local frame, the equilibrium velocity due to galac- 
tic differential rotation relative to the center of the box is given by Vq = —qQxy, where 
q = — d In Q/d In R\jif^ is the local shear rate and x = R — Rq and y = Ro{(f) — Qt) refer 
to radial and (rotating- frame) azimuthal coordinates, respectively. We adopt q = 1, cor- 
responding to a fiat rotation curve. The basic equations we solve are presented in Paper 
III. We use a modified version of the ZEUS code (Stone & Norman 1992a,b) with a hybrid 
Green function/FFT solver for the Poisson equation (Paper III) and applying a velocity de- 
composition method (Paper I) for less diffusive transport of hydrodynamic variables under 
strong shear conditions. For improved numerical performance in the low-density "coronal" 
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region, we employ the Alfven limiter algorithm of Miller &; Stone (2000) with the limiting 
speed of the displacement current, = Scg. We adopt periodic boundary conditions in the 
y- and 2;-directions and shearing-periodic conditions at the x-boundaries (HGB). From test 
runs with an inflow/outflow ;2-boundary, we have confirmed that evolution of the midplane 
portions of the disk is quite insensitive to the imposed vertical boundary conditions. 

The (x,y,z) dimensions of our computational box are {L^, Ly, L^) = (25,25,6) x H. 
Since an unmagnetized disk with Sq = 1 has H ^ 160 pc for the solar neighborhood value 
Q = 26kms^^kpc^^ (Paper 111), the total mass contained in the simulation domain is 
Mtot = 2.1 X lO^MG(Eo/13M0pc-2)(i//16O pc)^, which is about 20 times larger than the 
two-dimensional Jeans mass. 



The numerical resolution is 128^ zones for all models. We apply initial perturbations only 
to the velocity using spatially uncorrelated, random fluctuations multiplied by the back- 
ground density profile po{z). The amplitude of the perturbations is determined such that 
Max(|5v|/cs) = 0.25 at the disk midplane. We will use the orbital period ^orb = 27r/n = 
2.4 X lO^yrs (0/26 kms~^ kpc~^)~^ as the time unit in our presentation. 



To investigate MRI-driven turbulence and its interaction with self-gravity, we have per- 
formed 7 numerical simulations. The important parameters and the simulation outcomes 
for each model are listed in Table 1. Column (1) gives the label for each run. Columns (2) 
and (3) are the Toomre Q parameter and the initial field strength in terms of the plasma 
parameter at the disk midplane j3{0) (see eqs. [3] and [4]). Column (4) indicates whether 
self-gravity is included in the dynamical evolution. Note that self-gravity in the models with 
"yes" is slowly turned on only after turbulence becomes saturated, starting at t/torb = 4.1. 
Note also that "no" self-gravity models still include the initial total gravity that maintains 
the unperturbed equilibrium density profile. Column (5) gives the final simulation outcome: 
"unstable" indicates the formation of gravitationally bound condensations, while "stable" 
indicates fiuctuating density fields without any condensations. 

Columns (6) and (7) summarize the scaling relations among the time and space averaged 
stresses and magnetic energy density in the stable models: Rxy = {{pVxSvy)) and Mxy = 
{{—BxBy/An)) are the Reynolds and Maxwell stresses, respectively; T^y — Rxy + Mxy is the 
total stress; Eb = ((-B^/Stt)) is the total magnetic energy density. Here, the double bracket 
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(()) denotes a time- and space-average. Finally, columns (8)-(10) give the time-averaged 
values of density- weighted velocity dispersions {pSvfY^"^ / {p)^/"^ for the stable models, where 
the single bracket ( ) denotes a space average (with the subscript i = x, y, or z). For columns 
(6)-(10), time averages are taken typically over 2 orbits after MRI saturation. 

3.1. Non-Collapsing Models 

In this subsection we describe the evolution of non-self-gravitating and weakly self- 
gravitating models that do not form bound condensations. Of previous work, nonstratified 
simulations of HGB with vertical initial magnetic fields, or stratified disks of Stone et al. 
(1996, hereafter SHGB) with sinusoidally varying fields, are similar to our models. Both 
HGB and SHGB differ from the present models in that they used a Keplerian rotation 
profile appropriate to accretion disks. We shall compare the results of our simulations with 
those of HGB and SHGB. 



3.1.1. MRI-Driven Turbulence 

We begin with non-self-gravitating model A. Figure 1 plots evolutionary histories of the 
maximum density (pmax, solid lines) and density dispersion ((5p^)^/^, dotted lines) of model 
A together with those of the other models. Snapshots of the perturbed density at four time 
epochs of model A are shown in Figure 2. 

The initial equilibrium of model A is a vertically stratified, shearing disk with q = 1 and 
threaded by uniform Bz with /3(0) = 100. During the initial relaxation stage (t/torb ^ 2), 
small-scale modes have low-amplitude (a few percent) density fluctuations, as shown in the 
perturbed density plot (Fig. la). As the system finds the most unstable MRI mode, however, 
density fluctuations begin to grow rapidly. The dominance of the most-unstable axisymmet- 
ric mode (Balbus & Hawley 1991) is evident in the character of the density fluctuations of 
Fig. 2a, b, which are correlated with the "channel flow" in the velocity field (Hawley & Balbus 
1991; see below). Linear theory predicts that for the parameters in model A, the maximum 
growth rate of the MRI in model A is jmax = f^/2, which occurs for vertical wavelength 

A„„ . 200 PC (M) ""^ (^) {^^r^J^) (6) 

(e.g., Balbus & Hawley 1998). Since Amax fits well within a density scale height of the disk 
(i.e., Ainax/2 < H), the MRI in model A develops similarly to the way it would in a uniform 
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medium.^ The predicted maximum growth rate, 7max, marked by the dashed hne in Figure 
la, is indeed in good agreement with the simulation results. 

Because the fastest growing mode in model A approximates an exact nonlinear solution 
of the full MHD equations in the incompressible limit (Goodman & Xu 1994), it can grow 
to very large amplitude while preserving its essential character. As the mode grows, the 
whole flow in model A becomes dominated by this "channel solution," in which fluid ele- 
ments in alternating horizontal layers that have larger (smaller) azimuthal velocity than the 
background tends to move radially outward (inward). This inward/outward flow drags the 
magnetic fields radially in dircctionally alternating layers, and the sheared azimuthal motion 
produces corresponding alternating layers of toroidal magnetic field. As the amplitude of 
va,<p becomes large, vertical magnetic pressure gradients compress the disk in alternating 
layers. The gas associated with the over-dense regions at \z\ < O.OSL^ in Figure 2b is stream- 
ing radially inward, while the layers surrounding it are moving outward. After the channel 
mode grows to very large amplitTidc. the nonaxisymmetric parasitic instability (Goodman 
& Xu 1994) develops to break up the flow into smaller-scale turbulence (Fig. 2c, d). Direct 
growth of nonaxisymmetric MRI may also contribute to small scale kinetic and magnetic 
turbulence (HGB). At t/torb ~ 4, the turbulence has achieved a saturated state with large 
fluctuation amplitudes in all physical variables. 

The evolution of averaged disk properties due to the MRI is illustrated in Figure 3, 
where we plot horizontally-averaged densities, plasma /?, magnetic and perturbed kinetic 
energy densities, and Maxwell and Reynolds stresses of model A at times t/torb = 0,3,4,5. 
We also list in Table 2 the volume and time averages of various physical quantities averaged 
over t/torb = 4 — 6. Magnetic fields are amphfied to attain levels /3 ~ 1 — 10 at MRI saturation 
(Fig. 36). The regions above \z\ — 2H have /3 < 0.5, consistent with the results of SHGB and 
Miller & Stone (2000) that MRI in a stratified disk produces a moderately-magnetized disk 
midplane surrounded by a magnetically-dominated corona. Overall, the magnetic energy 
density is dominated by the azimuthal component (-B^/Stt), and is larger than the perturbed 
kinetic energy density by a factor of 2 to 7, with the largest value at late times occurring 
about a scale height above the midplane (Fig. 3c). Although the kinetic energy density is not 
more than ~ 10% of the central thermal pressure and the velocity fluctuations are subsonic 
at the midplane, the magnetic fields are strong enough to drive supersonic motions in the 
low-density, high-altitude regions. 

Tables 1 and 2 suggest that in terms of Shakura & Sunyaev (1973) a parameter de- 



•^The fact that Xmax/^z ~ 17, where is the grid spacing in the z-direction, means that the fastest 
growing wavelength is numerically well resolved. 
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scribing the radial transport of angular momentum, the total dimensionless stress is a = 
Txy/Po{0) ~ 0.16, with the Reynolds stress contributing only ~ 10% of the total (see also 
Fig. 3d). This demonstrates that magnetic fields bear primary responsibility for mediating 
mass accretion in an MRI-unstable disk (Balbus & Hawley 1998). During the entire linear 
and nonlinear evolution, the averaged density profile in z is within 20% of the initial dis- 
tribution, although the amplified magnetic field thickens the disk slightly via its magnetic 
pressure support (Fig. 3a). 

In spite of the kinematic effect of the background shear, the density and magnetic fields 
arc dominated by structure at relatively large spatial scales. Figure 4 displays the x-y density 
structure as well as horizontal velocity and magnetic field vectors at the midplane of model 
A at t/torb = 5.6. Fourier power spectra of the density^ and the azimuthal component of the 
magnetic field averaged over t/torh = 5 — 6 in model A are plotted in Figure 5. Although 
the volume-averaged values of and By are negligibly small (see Table 2), a constant-z 
plane consists of a few magnetic domains in which and By are distributed coherently; 
across the domains, the horizontal magnetic field reverses direction. For example, the right 
panel of Figure 4 shows that the midplane is dominated by negativc-i^j^ and positive- the 
midplane mean values of B^ and By are -0.18 and 0.38 times (47rP(0))i/2 t/torh = 5.6, 
respectively. Given these large mean values of (B^) and (By) in any layer at constant z, the 
magnetic field has to bend and reverse in the vertical direction in order to yield the small net 
volume averages of Table 2 {first three rows, second column). This field geometry is indeed 
consistent with Figure 56, which shows that the maximum power of By is attained in the 
Tlx — — riy, riz — 1 mode; the riy — 0, — 1 — mode is also strong. Here, integers 
Tlx = kx/27rLx, Uy = ky/27iLy, and = kz/27iLz denote the wavenumber in each direction. 
The dominant mode for By at late times is thus still that associated with channel fiow. The 
By mode with next most power is coherent in the y (azimuthal) direction and has a single 
reversal in the z and x (radial) directions. 

Despite the fact that our model disk has lower shear and is more tightly vertically 
compressed than the models considered by HGB and SHGB, the characteristics of turbulence 
observed in the power spectra are similar. Among the similarities are (1) most of power 
is contained in modes with scales comparable to or larger than the disk scale height, (2) 
the turbulence is highly anisotropic, as the different slopes of the power spectra along the 
different axes implies, and (3) at well-resolved scales, the logarithmic slopes of the power 
spectra are similar to the —11/3 Kolmogorov spectrum (see HGB; SHGB). The kinetic-to- 
magnetic scaling ratios in our models are also similar to those in HGB. For instance, for 



^Because the background density profile is stratified, there is a boost to the low-fc^ terms of the density 
power spectrum that is not associated with turbulent perturbations. 
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a nonstratified model (Z4) with uniform magnetic fields, HGB found Rxy/M^y ~ 0.19 and 
Txy/Ee ~ 0.66, which are very similar to R^y/M^y 0.13 - 0.14 and T^y/EB ~ 0.64 - 0.65 
in models A and B of this paper (see Table 1). 

The level of the turbulent fluctuations is, of course, dependent on the field strength 
and geometry, and possibly on other disk parameters. Having uniform vertical magnetic 
fields, model A exhibits stronger growth of MRI than a zero-nct-i?^ model (IZl) of SHGB. 
Table 2 shows that model A experiences fluctuations in the density with (5p^)^^^/(p) ~ 
1.018 and in the total pressure with {6P^^^.)^^'^ / (Ptot) ~ 0.629, which arc about an order of 
magnitude larger than in model IZl of SHGB. Comparing to the uniform- i?^ model Z4 of 
HGB, however, model A shows smaller turbulent fluctuation amplitudes. This is because 
the density stratification of model A makes high- 2; fields effectively stronger (in terms of /?), 
resulting in local MRI wavelengths at high latitudes that are larger than the vertical size of 
box. As a consequence, the resulting a in model A (which averages over both high-stress 
and low-stress volumes) is about a half oi a ^ 0.31 in model Z4 of HGB. 

In practice, the vertical magnetic field will always be correlated in sign on some largest 
radial scale, and the extent of this correlated region - as well as the magnitude of the vertical 
magnetic flux over this scale - will affect the amplitudes of saturated-state turbulence. Even 
over short evolutionary times compared to the galactic lifetime, our models show growth 
in the = 0, nj, = 0, n^; = 4 component of the Bz power spectrum, corresponding to 
correlated vertical magnetic flcld over radial scales comparable to the thickness of the disk. 
Although realistic values for the correlation scale and magnitude of arc not known at 
present, this evidence of "dynamo" activity helps to motivate our incorporation of vertical 
flux in the initial conditions. While there are not yet predictions from simulations of the 
level of vertical magnetic flelds that would develop over long timescales, one may speculate 
that vertical fleld would build up to values similar to those we adopt, since for much weaker 
vertical flelds the fastest-growing MRI modes are at very small scales, and for much stronger 
vertical flelds the MRI is stabihzed (see also §4). 

Radio observers of external face-on galaxies report line widths of order ~ 6 — lOkms"^ 
in extended regions of H I disks where no active star formation takes place (Dickey et al. 
1990; van Zee & Bryant 1999). Scllwood and Balbus (1999) proposed that in the absence of 
stellar energy inputs, the MRI may be important in producing large-scale random motions in 
these extended H I disks. ^ A direct way to assess this specific proposal, and more generally 
to investigate the role of the MRI in driving galactic turbulence, is to measure the turbulent 



^An alternative explanation for turbulence in extended disks without stellar energy feedback and in NGC 
2915 attributes it to self-gravitation and bar driving (Wada et al. 2002); see also §4.1. 
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velocity dispersions from numerical simulations. Although the current models are idealized 
in ways that may significantly affect the results, it is useful to quantify the turbulence 
level for future comparisons. As Table 1 lists, the density-weighted one- dimensional velocity 
dispersions in the non-self-gravitating model A are found to be a^j. = 0.45Cs, ay = 0.43Cs, and 
(T2 = 0.23cs. In dimensional units, these correspond only to ~ 1.6 — 3.2 kms~^, much smaller 
than the observed line widths cited above (which however are not purely turbulent - see §4.1). 
The velocity dispersions are smallest in the vertical direction, indicating that buoyancy does 
not appreciably enhance turbulence. Model Z4 of HGB reports (Jx = (Ty = OAScg and 
az — 0.22cs, so that the density stratification has httle effect on the velocity dispersion. 

Finally, we compare the evolution of our disk models with Miller & Stone (2000), who 
performed stratified disk models with q = 3/2. In their simulation with a pure field, 
strong axisymmetric channel fiows quickly emerge and rise buoyantly, disrupting the initial 
disk structure. The system becomes magnetically dominated everywhere at the end of the 
run. As Figure 3 shows, however, our model A keeps its initial disk structure fairly intact 
even when the disk is highly turbulent. High-density parts of the disk in model A have 
subthermal magnetic fields at saturation rather than becoming magnetically dominated. In 
part, these differences may be due to the differences between the circumstellar-disk tidal 
gravity law = —fl'^z and the stronger vertical gravity provided by the gaseous and stellar 
disks in a galaxy. Within a scale height of the midplane, the galactic gravity is a factor 
~ 2 — 4 larger for the present models compared to — —Q^z. This compresses the disk 
more tightly in the vertical direction, making the buoyant rise associated with the channel 
solution less efficient. With an initial field strength having /3(0) = 12.5 and a shear parameter 
q — 3/2, MRI-amplified channel-flow magnetic fields in Miller & Stone (2000) are also much 
stronger than those in the present models, contributing to the destruction of the initial disk 
structure. 



3.1.2. Effect of Self- Gravity on Turbulence 

We now turn to models F and G, which have self- gravity included in their dynamical 
evolution but have surface densities below the (nonlinear) gravitational instability threshold. 
Thus, self-gravity in models F and G, with Q = 1.7 and 2.0, respectively, is not strong enough 
to engender formation of bound condensations, although it certainly enhances the fluctuation 
amplitudes of physical variables. Table 3 lists the volume and time averaged quantities (over 
t/torh = 5.0 — 6.5) in model F. The fluctuations in the perturbed density and total pressure, 
with {Sp'^y^^/{p) ^ 1.021 and {SP^^^y^"^ / (Ptot) ~ 0.646, are quite similar in amplitudes 
to those in model A, but magnetic and perturbed kinetic energy densities in model F are 
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increased by about 60%. Note, however, that the ratios Rxy/M^y and Txy/EB are almost 
independent of the presence of self-gravity (Table 1), indicating that self-gravity, as well as 
the disk compression and shear parameter, does not change the angular-momentum-transport 
characteristics of MRI-driven turbulence significantly. 

To gain some idea on the origin of the extra turbulent energy in the (weakly) self- 
gravitating model F, wc have analyzed for comparison the data of the hydrodynamic (i.e. 
unmagnetizcd) model C (with Q = 0.7, Sq = 25, and f3 = oo) of Paper 111. The density- 
weighted velocity dispersions, averaged over t/torb = 1-2 — 2.2, are found to be = 0.12cs, 
a-y = O.OScg, and cr^ = O.OScg, which happen to be similar to Aa^ = O.lScg, Aay = 0.07cs, and 
Aaz = 0.07cs, the differences in the velocity dispersions between the MRl-unstablc models F 
and A. Note that the turbulence in model C of Paper III is driven purely by hydrodynamic 
swing amplification (perhaps augmented by buoyancy effects). Although it is not obvious 
how turbulent energies of different physical origins in general ought to combine to reach 
a total turbulent amplitude, these results suggest that motions driven by self-gravity in a 
single-phase medium can contribute only ^ 20% of the amplitude of the velocity fluctuations 
that MRI can drive. Therefore, self-gravity alone^ is unlikely to be a significant source of 
turbulence in galactic disks with Q ^ 1.7. Disks with smaller Q become gravitationally 
unstable, as we shall show below. 



3.2. Self-Gravitating Cloud Formation 

3.2.1. Q Threshold for Runaway and Bound Cloud Masses 

We now consider the formation of self-gravitating clouds in turbulent galactic gas disk 
models. As mentioned before, we slowly turn on self-gravity at t/torb = 4.1, after MRI 
saturation. This slow inclusion of self-gravity prevents an abrupt response to the gravi- 
tational potential and enables us to explore the nonlinear interaction of self-gravity with 
fully-developed turbulence. 

As Figure 1 shows, models C-E with Q < 1.5 achieve high enough density to pro- 
duce self- gravitating runaway, while other models with Q > 1.7 remain stable (with large 
fluctuations in density). From careful examinations of density power spectra at successive 
times, we identify the mechanism of self-gravitating cloud formation in models C-E with 
swing ampliflcation of high-amplitude density perturbations. As discussed in the previous 



^Scc however the discussion in §4.1 of gravitational driving of turbulence when thermodynamic evolution 
can maintain a very low filling-factor medium. 



-14- 



subsection, the power spectra of MRI-driven turbulence peak at small wavenumbers, so that 
the corresponding large spatial scale density perturbations can be easily swing-amplified if 
self-gravity is sufficient. Figure 6 shows surface density snapshots at t/torb=4.5, 5.0, and 5.5 
of model C (with Q = 1.0). In the left frame, the maximum power in the surface density is 
in a trailing rix = 4:,ny = 1 mode, but there is also considerable power in a leading mode 
with Hx = —I and Uy = 1. The swing amplifier causes the leading mode to accumulate 
more mass as it shears into trailing configuration. The net result is the formation of high 
density filaments (middle frame), which collide with each other and then fragment to form 
three bound condensations (right frame). At their initial formation, the clumps have mass 
M Mj_2D ~ 10^ Mq each (see eq. [5]); they subsequently accrete additional material from 
their surroundings. At the end of the simulation (t/torb=5.5), the mass contained in the 
three clumps is ~ 14% of the total mass, corresponding to each clump of mass ~ 2.2Mj 2d- 

Because of its weaker self-gravity, model D with Q — 1.5 and /3(0) = 100 requires two 
major swing amplification events before its eventual gravitational runaway occurs. In partic- 
ular, intermediate-density filaments (from the first swing amplification) have Pmax ~ 4.6po(0) 
at t/torb ~ 5, but these fail to initiate gravitational collapse. Although most of power at 
this point is in the form of trailing waves, nonlinear interactions of sheared large-amplitude 
wavelets are vigorous, especially given the stresses from highly turbulent magnetic fields, 
which ends up supplying fresh small- |A;a;| modes (Paper 1). Another phase of swing amplifi- 
cation follows, this time yielding sufficient enhancement of density to trigger gravitational 
runaway. Model D forms a bound clump with M Mj_2d at the end of the run. The 
evolution of model E with Q — 1.5 and /3(0) = 400 is quite similar to that of model D, but 
the lower-amplitude saturated state of MRI when the mean magnetic field is weaker causes 
the formation of a bound cloud to take longer. 

The self-gravitating clouds that form in our simulations turn out all to be magnetically 
supercritical. For instance, the mass-to-fiux ratio of the cloud marked by a square in the left 
frame of Figure 6 is M/^b ^ 1.3G^^^^ at t/torb = 5.1 when it first appears, and increases 
to ~ 3.6 G-^/^ at t/torh = 5.5, weU above the critical value of ~ 0.13 G^^/^ (e.g., Shu 1992). 
Here, $5 denotes the magnetic fiux that passes through the cloud and G is the gravitational 
constant. The other bound clouds also have mass-to-fiux ratios in the similar range, that 
is, M/^B ~ (1 — 4) with higher values corresponding to more-evolved, more-massive 

clouds. Note that the lower limit of this range is remarkably similar to M/^b ~ 0.96 G"^/^ 
of a square region projected in the x — y plane that contains one Jeans mass in the saturated- 
MRI disk with the mean By ~ 2/iG. This suggests that the bound clouds initially form by 
collecting material relatively isotropically in a way that preserves the magnetic fiux fairly 
well. As the clouds become sufficiently self-gravitating, they begin to collapse and accrete 
more gas in all three dimensions, perhaps with accretion along the mean field direction 
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enhancing the mass-to-flux ratio at later time. Although hard to quantify with our present 
limited resolution, magnetic reconnection occurring inside the clouds may also contribute to 
high mass-to-flux ratios. 

As Figure 1 clearly illustrates, the threshold for nonhnear gravitational instability in our 
MRI-driven turbulent galactic disk models is in the range Q — [1.5, 1.7]. This threshold Qth 

is higher than found for either completely unmagnetized {Qt\i < 1) or thermal-equipartition- 
magnetized (Qth ~ 1) 3D disk models in Paper III. When bound clouds formed in those 
models, they had masses of two to several Mj^2Di similar to the present results. We note 
that in our previous simulations without large-amplitude turbulence, any model with density 
fluctuations as large as those in model F (see Table 3) ended in gravitational runaway. Here, 
however, the same turbulence that creates overdense regions can also destroy them, such that 
model F is nonhnearly gravitationally stable. The MRI-driven turbulent disks have higher 
Qth than the models of Paper III partly because the turbulence produces large-amplitude, 
large-scale perturbations in the surface density (see Tables 2,3) that reduce the effective 
value of Q by 20 — 40%. In addition, torques from small-scale magnetic fields can drive 
radial angular momentum and mass transport, much as occurs from the MRI in accretion 
disk models. Since vorticity conservation of fiuid elements is the chief barrier to gravitational 
instability at large enough scales in rotating disks, radial transfer of angular momentum 
within/across an overdense region can contribute to enabling runaway condensation (cf. 
Gammie 1996). 



3.2.2. Angular Momenta of Condensations 

A longstanding astronomical issue of great interest is how the total angular momentum 
of a self-gravitating cloud evolves as it condenses out of the diffuse ISM. The prevailing idea 
is that magnetic fields linking a cloud with the surrounding medium exert back-torques as 
it spins up during contraction, leading to loss of angular momentum (e.g., Gillis, Mestel, & 
Paris 1974; Mouschovias & Paleologou 1979). Most previous work on this "magnetic braking" 
process has focused on the problem of an ideal cloud anchored to the ambient medium with 
a prescribed magnetic field distribution, solving to obtain the temporal evolution of the 
angular velocity. Using our three-dimensional simulation data from both the present models 
and those from Paper III, we will demonstrate directly in this subsection that a cloud that 
condenses out of a rotating, magnetized disk can shed significant angular momentum via 
magnetic braking during the cloud-forming phase. 

To describe the rotational properties of a self-gravitating cloud, we consider the volume 
bounded by an isodensity surface p = Pcrit- We define the total angular momentum and the 
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total mass of a volume with p > pcrit as 



(7) 



P>Pcrit 



A^(Pcrit) = / Pd^x, 



(8) 



P>Pcrit 



respectively, where the coordinates (x, z) and the inplane velocity (f^, f^^) are measured 
relative to the center of the condensation. One may then calculate the mean specific angular 
momentum = Lzj hA as a function of Pcrit- Because model disks from which bound clouds 
originate may have different specific angular momenta initially, it is useful to compare 
of a cloud with Jgai, the total specific angular momentum that the cloud would have if it 
did not lose any angular momentum during its formation. For a region of uniform surface 
density, it can be shown that 



where A is the area of the region in the plane of the disk (cf. Mestel 1966). Thus, in the case 
of a fiat rotation curve = 1, Jgai = ^L\lVl for a square patch of side Lq in the unperturbed 
disk. For purposes of comparison, we shall take Lq such that the total mass is the same as 
that within the p — pcrit surface for the bound cloud that forms. 

Exemplary profiles of JziPcrit) / Jgai ^oi three bound clouds at a few selected times are 
shown in Figure 7. In the figure, "MRI" refers to a cloud that forms within the region 
marked by a square in the Figure 6 snapshot of model C {Q — 1.0, sq = 1) of this paper, 
while the curves with "Parker" and "Swing" correspond to clouds formed in magnetized 
model D ( Q = 0.7, Sq = 25,/9 = 1) and unmagnetized model C {Q = 0.7, Sq = 25,f3 = oo) 
of Paper III, respectively. It is apparent that magnetized clouds ("MRI" and "Parker") 
lower their specific angular momentum over time as they gravitationally contract, while Jz 
of an unmagnetized cloud ("Swing") fiuctuates around ~ O.SJgai- This implies, consistent 
with expectations, that magnetic fields play a key role in removing angular momentum. 
Although the resolution of the present simulations is not sufficient to follow this process over 
long periods, the characteristic angular momentum loss time during the interval measured is 
Jz/\Jz\ ~ (1 — 2)Q~^ for the magnetized models. Since bound clouds typically form within 
^cou ~ (1 — 3)J1~^ (see Fig. 1), the similarity of timescales suggests that magnetic braking is 
quite efficient, possibly enough to explain the very low observed specific angular momenta 
of GMCs (e.g., Bhtz 1993) compared to their corresponding Jgai values. 

To help visualize the magnetic braking process that occurs in our simulations, we show 
in Figure 8 a volumetric rendering of an isodensity surface and magnetic field lines, as well as 




(9) 
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the density and velocity field at the midplane, for a dense clump in model C at t/torh = 5.5. 
The midplane velocity vectors show that the cloud is counter-rotating with respect to the 
the sense of the background shear, with the principal rotation axis parallel to the 2;-direction 
(this sense of spin is prograde with respect to galactic rotation) . 

The magnetic field lines show rather complex behavior. Some field hues (black) at 
mid-latitude regions loop back to the cloud, while other field lines anchor the cloud to the 
ambient medium. Field lines drawn in green start running from the right edge of the box, 
swirl into the cloud in the counterclockwise direction, and rccmcrgc to swirl in the clockwise 
direction out to the left edge of the box. These field lines exert net torques on the cloud 
surface. On the other hand, field lines drawn in blue that touch the ceiling of the box 
are relatively straight. The fact that the horizontal field lines are more twisted than the 
vertical field lines suggests that magnetic fields perpendicular to the spin axis provide most 
of the braking. Although the periodic ^-boundary conditions we have adopted tend to reduce 
vertical torsion, an important physical consideration is that there is larger inertia of ambient 
gas in the horizontal direction (e.g., Mouschovias & Paleologou 1979), tending to promote 
large horizontal torsion. 

Finally, we remark that limited resolution in the current simulations may affect the 
measured estimate of angular momentum loss. With higher-resolution simulations (employ- 
ing adaptive or nested meshes), it will be possible to make a more quantitative assessment 
of how effective magnetic braking is in removing a cloud's initial angular momentum. 

4. Summary and Discussion 

Understanding the origins of GMCs is key to characterizing star formation on a galactic 
scale, because the rate and mode of star formation are linked to the rate of cloud formation 
and the physical properties of GMCs at birth. Many different galactic-scale processes may 
affect GMC formation. In this paper, we use 3D MHD simulations to initiate investigation of 
MRI-driven turbulence in disk galaxies and its role in prompting cloud-forming gravitational 
instabilities. 

Our numerical models represent local portions of shearing and vertically stratified gas 
disks with initial density profiles determined by the balance between thermal pressure gra- 
dients and gaseous self- and external stellar gravity (see §2). The background azimuthal 
velocity is set to have a local shear rate of g = —dhiVL/dhiR = 1, corresponding to a fiat 
rotation curve. Magnetic fields are initially assumed to be purely vertical and uniform, with 
strength characterized by /3(0) = 100 or 400 (see eq. [3]). For simpficity, we have adopted 
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an isothermal equation of state with sound speed Cg = 7kms~^; the corresponding temper- 
ature is characteristic of the warm ISM. We evolve disk models of varying surface density 
parameterized by the Toomre stability parameter Q (see eq. [4]); we consider models with 
Q ~ 1—2. We explore the nonlinear saturation of MRI and the growth and evolution of 
self-gravitating structures. In particular, we analyze statistical properties of the turbulence, 
assess the threshold for gravitational instability, and measure the masses and angular mo- 
menta of the bound clouds that form. Where appropriate, for comparison we also analyze 
the data from some models from Paper III. 



4.1. Properties of Turbulence 

In §3.1 wc investigate the characteristics of turbulence initiated and sustained by MRI 
in non- or weakly self-gravitating models. Despite the differences in the shear parameter, 
degree of vertical disk compression, and the presence of self-gravity, the statistical properties 
of the turbulence in our stratified disks are remarkably similar to those in the nonstratified 
"accretion disk" simulations of HGB. Our models yield saturated-state Shakura & Sunyaev 
dimensionless R—^p stress parameter a ~ 0.15 — 0.3, ratio of the Reynolds to Maxwell stresses 
Rxy/Mxy ~ 0.13, and ratio of the total stress to magnetic energy density T^y/Es ~ 0.64 (see 
Table 1). The shapes of the turbulent power spectra are also similar to those of HGB. With 
a level a ~ 0.2 (see also Hawley & Krolik 2002) characterizing angular momentum transport 
and Cg = 7kms~^, the MRI-driven gas accretion time t^cc ~ R'^fl/{acl) would exceed 10^^ 
yr at the solar circle. 

To quantify the level of turbulence that is driven by MRI in our models, we directly 
measure the one-dimensional velocity dispersions in each coordinate direction. We find 
density-weighted velocity dispersions in models with /3(0) = 100 are ax ^ <Jy ^ (0.4 — 0.6)cs 
and az ~ (0.2 — 0.3) Cg, with the mean plasma parameter /9 ~ 1 — 2 at saturation. With 
Cg = 7kms~^ and po(0) ~ 10~^^ g cm~^, this implies that MRI in the present models 
produces turbulent levels up to cTj; ~ (t^^ ~ 4kms~^ and ~ 2kms~^, and RMS magnetic 
field strength S ~ (1.8 — 2.3)//G. These magnetic field strengths are similar to estimates in 
the Milky May and in external galaxies (e.g.. Rand & Lyne 1994; Beck 2000, 2002). The 
velocity dispersions are significantly lower than the observed line widths of ~ 7kms^^ and 
~ llkms^^, respectively, for the cold and warm H I components in the Solar neighborhood 
(e.g., Heiles & Troland 2003). Since the thermal linewidth from cold gas is negligible and 
that from the warm gas (in the stable regime) is ~ 7kms^^, the warm and cold phases both 
have turbulent velocity dispersions ~ 7kms~^, significantly larger than the turbulence levels 
in the present models. For the extended disks of face-on galaxies, radio observations of H 1 
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emission linewidths imply (j^ ~ 6 — lOkms"^ (e.g., Dickey et al. 1990; van Zee & Bryant 
1999; Petric & Rupen 2001) for the atomic gas; in the absence of absorption measurements, 
it is not directly possible to distinguish the turbulent and thermal contributions (or the 
proportions of warm/cold gas).^ 

Because the level of the velocity fluctuations depends on the mean magnetic flux, one 
might argue that we could have obtained larger a's, had we adopted larger initial vertical 
5-fields. HGB found, for example, that the saturation magnetic energy density is linearly 
proportional to Amax and thus to the mean value B,{0), or P{0)-^/^. If the MRl in a single- 
phase medium is to produce up to the observational values, these scalings would imply 
that the vertical magnetic field should have j3{0) ^ 10 (assuming that the magnetic and 
perturbed kinetic energy densities at saturation are also linearly proportional). Equation (6) 
would then require Amax ~ 600 pc for the vertical wavelength of MRI; if this is larger than the 
thickness of the disk, the saturation energy density would have to be reduced accordingly. 
Even if larger mean vertical fields were able to produce larger velocity dispersions, the implied 
total magnetic energy densities would exceed observed values in disk galaxies. Given these 
constraints, we conclude that MRI acting in a single-phase medium could not generate the 
level of turbulence implied by H I line width measurements. 

It is unknown - but of great interest - whether MRI processes acting in a much less 
uniform, multi-phase disk could yield appreciably larger velocity dispersions. One might 
expect, for example, a larger velocity dispersion for a given MRI driving rate if larger effective 
mean free paths for dissipative interactions in a cloudy medium (compared to a single-phase 
medium) reduces the turbulent decay rate. Even if the MRI angular momentum transport 
rate - proportional in steady state to the energy dissipation rate per unit mass a{P/p)qfl - 
increased by a factor 4 — 9 for a factor 2 — 3 increase in (Jturb compared to the above, the 
accretion time (oc 1/a) is still comparable to the Hubble time. This timescale for inflow 
is long enough to allay the historical concern of overrapid accretion if turbulence taps the 
orbital energy of the ISM in a galaxy (cf. Spitzer 1968). In addition, the implied volume 
heating rate would also still be well within the limits set by C II emission from the diffuse ISM 
(Wolfire et al. 2003). Numerical simulations will be required to decide whether MRI-driven 
turbulence in a multi-phase medium is competitive with the conventional turbulent energy 



^Interestingly, if cr^^^ = cr^uj-b + /warmC^ warm where /warm is the mass fraction in the warm phase (and 
assuming CTturb is the same for all gas), then values of Uobs lower than Cs,warm ~ 7kms~^ suggest that (1) 
/warm ^ 1, i.e. at least some cold gas is present (consistent with theoretical expectations, cf. Wolfire et al. 
2003), and (2) CTturb is lower than the ambient medium's sound speed Cg. warm- For example, ifcrturb ~ 2kms~"'^ 
as found for cr^ in our models, then with Cg^warm ~ 7kms^^, dobs = 6kms~^ would imply /warm = 0.65 and 
/cold = 0.35. 
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source for the ISM, supernova shocks (e.g., Cox & Smith 1974; McKee & Ostriker 1977). 
Initial work towards this goal is already underway (Piontek & Ostriker 2003). Realistic 
theoretical assessment of the contributing turbulent driving processes may ultimately require 
simultaneous modeling, since the dynamical effects are not necessarily independent. 

Our three-dimensional simulations show that the enhancement in the amplitudes of ve- 
locity fluctuations by the presence of self-gravity is so small (less than 30%) that self-gravity 
alone is unlikely to be a major source of galactic turbulence for a single-phase disk. Recently, 
Wada & Norman (1999) and Wada et al. (2002) using two-dimensional hydrodynamic simula- 
tions argued that self-gravity can extract sufficient turbulent energy from large-scale galactic 
rotation to maintain the observed level of turbulence in the ISM. Their conclusion is based 
on the development of a turbulent velocity field with (Jturb ~ lOkms"^ in thermally and 
gravitationally unstable flows. At the point a quasi-steady state is reached in the Wada et 
al. models, most of the gas ( ^ 80%) is contained in cold {T < lOOK) clouds occupying < 1% 
of the disk area. Thus, the cloudy medium resembles a disk of collisionless particles which 
has kinetically heated itself through large-scale gravitational instabilities until the typical 
dispersion in random velocities ~ lOkms"^ is sufficient to render Qefr > 1- These sorts 
of models are likely appropriate for the ISM in galactic center regions (often dominated by 
molecular gas) rather than for the ISM in the main and outer disk, where surface densities are 
an order of magnitude lower than those considered by Wada et al., and where observations 
indicate that a much smaller proportion of the gas is in cold, dense clouds. 

In an application to NGC 2915, Wada et al. (2002) also presented a 2D model that 
extends up to i? = 15 kpc. They found that rapid cooling of gas makes even very low 
surface-density disks highly susceptible to thermal and gravitational instability that produces 
a cloudy medium with (Jturb ~ 2 — 4kms~^. Combined together with our finding that turbu- 
lence driven by self-gravity is not significant as long as the disks are weakly self-gravitating 
(with Q ^ 1.7) and remain isothermal, this suggests that the dynamical importance of self- 
gravity may depend strongly on the effective filling factor of the gas. It is therefore essential 
to run 3D models with a range of properties to determine how the amplitudes of turbulence 
from self-gravitational driving scale with the relative proportions of cold, dense, and warm 
(or hot) diffuse gas. These proportions of gas in different phases cannot be tuned directly 
since turbulence itself can affect the ratio of dense to diffuse gas, via coUisional shock heating 
and other related processes. 
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4.2. Bound Cloud Formation 

In §3.2, we study directly the nonlinear interaction of MRI-driven turbulence with self- 
gravity, showing that Jeans-mass-scale bound condensations form provided that the mean 
gas surface density is sufficiently high. In terms of the Toomre parameter, the nonlinear 
instability condition is Q < Qth ~ 1-6. The route to cloud formation is swing amplification 
of nonlinear density fluctuations 2.3 — 4.6) over large spatial scales (~ kpc) that 

arise as a consequence of the MRI.^ Compared to our previous 3D models of unmagnetized 
or strongly-magnetized (/9 = 1) disks which yielded Qth ~ 1 (Paper III), the MRI-unstable 
disks have larger Qth- This higher Qth may partly be because self-gravity is initially strong in 
the large-amplitude density fluctuations driven by turbulence, enhancing the effectiveness of 
the swing amplifler; the higher Qth may also owe in part to the ability of small-scale magnetic 
flelds to transfer angular momentum between over- and under-dense regions locally. 

The Q threshold we flnd is similar to the empirical result Qth ~ 1.4 obtained by Martin 
& Kennicutt (2001) based on azimuthally- averaged gas surface densities in a large sample 
of galaxies. Although intriguing, we caution that this close coincidence may be particular 
to the details of the current models (e.g. the net poloidal magnetic flux, the phase state 
of the gaseous medium). Nevertheless, we consider it a significant success that theory and 
observation have arrived at essentially the same simple criterion for the onset of active star 
formation in galactic disks. 

Gravitationally bound clouds that form in our simulations are all magnetically super- 
critical with the mass-to-flux ratio of M/^b ~ (1 — 4)G~^/^ and have typical masses of a 
few 10^ M0 each, corresponding roughly to the two-dimensional Jeans mass at the mean 
surface density. These clouds have masses consistent with the largest atomic or molecular 
clouds found in disk galaxies (e.g., Elmegreen & Elmegreen 1983; Vogel et al. 1988; Rand & 
Kulkarni 1990; Sakamoto et al. 1999). The molecular components of the largest Milky Way 
GMCs are nearly an order of magnitude lower in mass. The added contribution from H I, 
which increases total masses of Milky Way GMCs by at least a factor of two (Blitz 1993), 
reaching ^ 10^ Mq in some cases (e.g., Elmegreen & Elmegreen 1987), accounts for some of 
the disparity between our present results and observed molecular clouds. In addition, the 
disparity may owe in part to expected differences between cloud masses in systems with and 
without strong spiral structure, as discussed below. 



°In principle, turbulence driven by other processes - such as supernovae - would also be able to produce 

the high-amplitude density fluctuations needed to seed swing at Q-values as large as 1.5. However, it is 
not yet clear whether other mechanisms can produce density enhancements over the required ^ kpc spatial 
scales. 
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Table 4 summarizes the model prescription/parameters and main results for bound cloud 
formation from Papers I-III and the present work. As Table 4 indicates, our simulations so far 
have considered both disk models without features in the background gravitational potential 
(Paper I in 2D; Paper III and this work in 3D), and disk models that include a background 
(stellar) spiral potential (Paper II in 2D). All of our models to date are local. 

Although restricted to the 2D razor-thin limit. Paper II demonstrated that when spiral 
structure is present, gravitational instabilities (in particular, the MJI) form clouds preferen- 
tially in spiral arm regions, consistent with the observations that the warmest, most massive 
GMCs in the Milky Way arc strongly correlated with spiral structure (e.g., Solomon ct al. 
1985; Dame et al. 1986; Solomon & Rivolo 1989). While the mean gaseous surface densities 
for models that formed clouds in Paper I (without spiral structure) and Paper II (with spi- 
ral structure) were similar, the resulting cloud masses in Paper II were smaller by a factor 
~ 5 — 7. This difference owes to the overall compression of diffuse gas in spiral arms, which 
reduces the Jeans mass oc E~^; the mass collection zone is also limited by the width of the 
spiral arm itself. 

For a given gaseous surface density, cloud masses are also affected by the thickness of 
the disk, since increasing the thickness dilutes self-gravity, typically increasing the in-plane 
Jeans length by a factor ~ 2. With weaker self-gravity, a larger total area/mass is needed to 
prompt instability. Thus, the 3D models of Paper III produced larger cloud masses at given 
surface density (i.e. Q value) than would have been predicted on the basis of the 2D models 
of Paper I (although since instability was present only for smaller Q in 3D than in 2D, the 
resulting cloud masses in Papers I, III were similar). 

Real galactic disks of course have finite thickness and arc turbulent, and in this sense 
the current models supersede the "swing" models of Papers I and III. However, real galaxies 
also in general contain stellar spiral structure, which Paper II implemented only in 2D. 
Finite thickness effects are expected to increase the cloud spacings and masses that form by 
MJI in spiral arms compared to the results of Paper II. Although this could imply up to a 
factor ~ 4 increase in cloud masses at a given surface density, turbulence could limit this 
increase somewhat. Three-dimensional simulations with explicit modeling of spiral potential 
perturbations and vertical density stratification arc required to address this and related 
questions directly, although our net expectations from effects currently studied suggest total 
cloud masses in the neighborhood of 10^ M©. 

Because observed GMCs have very small specific angular momenta (Blitz 1993) com- 
pared to galactic values (see eq. [9]) on comparable mass scales, it is of much interest to 
determine when and how this angular momentum is lost. Using our models, we explicitly 
demonstrate that a self-gravitating, contracting, magnetized cloud loses specific angular mo- 
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mentum via surface torques imposed by the large-scale galactic magnetic fields that thread 
it. This magnetic braking is predominantly from horizontal S-fields (perpendicular to the 
spin axis of the cloud) that connect the cloud to the dense portion of the galactic disk. Our 
preliminary estimates suggest that the rate of angular momentum loss is sufficient to explain 
observations of low GMC spins, although higher-resolution simulations will be needed to 
confirm this result. 

Finally, wc remark on the mass spectrum of Galactic molecular clouds, generally de- 
scribed from CO observations as a power-law dN/dM oc M~^'^ below a cutoff of a few 
xlO^M0 (e.g., Solomon et al. 1987) (with H I envelopes increasing individual cloud masses 
by factors of 2 — 10). The gravitational instability scenario of this work and Papers I-III 
explains the upper cutoff in terms of the Jeans mass, as noted above (see also e.g. Elmegreen 
1979, 1994). Below this cutoff, the mass spectrum may mainly reflect the effects of turbulent 
dynamics within massive (Jeans-scale) clouds as they form and subsequently disperse. Re- 
cent analysis of simulated models of turbulent GMCs has shown that a power-law spectrum 
of moderate-density clumps develops (Ostriker 2002); the index measured is the same as 
the observed Galactic cloud spectrum.^ If the processes that disperse a GMC preserve the 
mass spectrum of former clumps as isolated clouds, and/or if observational identification 
methods decompose some large GMCs into their sub-parts, then similar clump and cloud 
mass spectra are a natural consequence of "top-down" GMC formation. 

It is a pleasure to acknowledge valuable discussions with S. Balbus, J. Dickey, B. 
Elmegreen, L. Mundy, R. Narayan, F. Shu, and S. Vogel. We also are grateful to the 
referee, K. Wada, for a insightful and constructive report. This work was supported in part 
by NASA grants NAG 5-9167 and 5-10780 and NSF grants AST 0205972 and AST 0307433. 
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Table 1. Parameters of Three- Dimensional Simulations 



Model 


Q 


m 


Gravity 


Results 






Oxlc^^ 






(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


A 


1.5 


100 


no 




0.13 


0.64 


0.45 


0.43 


0.23 


B 


1.5 


400 


no 




0.14 


0.65 


0.21 


0.21 


0.13 


C 


1.0 


100 


yes 


unstable 












D 


1.5 


100 


yes 


unstable 












E 


1.5 


400 


yes 


unstable 












F 


1.7 


100 


yes 


stable 


0.14 


0.62 


0.58 


0.50 


0.30 


G 


2.0 


100 


yes 


stable 


0.13 


0.64 


0.58 


0.50 


0.32 



'^Rxy and Mxy are the x-y components of the volume- and time- averaged Reynolds and 
Maxwell stresses, respectively. 

^Txy = Rxy + Mxyi and Eb is the volume- and time- averaged value of the magnetic energy 
density. 

^(Tj;, (Ty, Oz are the volume- and time-average of density- weighted velocity dispersions in the 
x-, 7/-, 2;-directions, respectively. 
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Table 2. Volume- and Time- Averaged Quantities in Model A 



Quantity 


// -f w 

((/)) 


/ /X f2\\l/2 


mm / 


max / 


S^/(47rPo(0))V2 


2.39x10"^ 


0.326 


-1.362 


1.317 


y/ \ ^\ / / 


-1.30 X 10"^ 


0.576 


-1.550 


1.467 


B /f47rPnfO))^/^ 


0.1000 


0.174 


—0.898 


1.189 


Bl/8nPo{0) 


0.0532 


0.072 


0.000 


0.957 


p2/87rPo(0) 


0.1808 


0.147 


0.000 


1.218 


p2/87rPo(0) 


0.0196 


0.035 


0.000 


0.691 


-P,P^/47rPo(0) 


0.1435 


0.168 


-0.880 


1.323 


-P,P,/47rPo(0) 


-0.0012 


0.075 


-0.850 


0.838 


-P^P,/47rPo(0) 


0.0031 


0.110 


-0.885 


0.913 


pvl/2Po{0) 


0.0334 


0.062 


0.000 


1.557 




0.0293 


0.065 


0.000 


1.687 




0.0116 


0.023 


0.000 


0.718 


pVjVy/Po{0) 


0.0180 


0.088 


-1.011 


1.760 


pv^vjPo{0) 


0.0005 


0.050 


-1.250 


0.883 


pSvyvjPo{0) 


0.0008 


0.050 


-1.050 


1.090 


Ptot/Po(0) 


0.5863 


0.396 


0.029 


2.478 


p/Po(0) 


0.3320 


0.338 


0.000 


2.184 


(P-Po)/Po(0) 


0.0000 


0.184 


-0.851 


1.233 


s/s„ 


1.0000 


0.201 


0.471 


1.779 
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Table 3. Volume- and Time- Averaged Quantities in Model F 



Quantity 


((/)) 


((5/2))l/2 


min / 


max / 


S2/87rPo(0) 


0.0767 


0.108 


0.000 


1.879 


S2/87rPo(0) 


0.3138 


0.215 


0.000 


2.076 


S,V87rPo(0) 


0.0230 


0.045 


0.000 


0.966 


-B,By/47iPo{0) 


0.2252 


0.249 


-1.378 


2.538 




-0.0008 


0.104 


-1.351 


1.412 


-ByBjAirPoiO) 


0.0030 


0.164 


-1.468 


1.483 




0.0566 


0.117 


0.000 


3.088 


p6v'y/2Po{0) 


0.0412 


0.111 


0.000 


3.398 


Pv!/2Po{0) 


0.0167 


0.033 


0.000 


0.888 


pVjVy/Po{0) 


0.0310 


0.156 


-2.307 


3.544 


pv^v^/Po{0) 


0.0023 


0.082 


-1.491 


1.928 


pSvyvjPo{0) 


0.0013 


0.074 


-1.669 


1.790 


Ptot/Po(0) 


0.7461 


0.482 


0.046 


5.474 


P/Po(0) 


0.3320 


0.399 


0.000 


5.068 


(p-Po)/po(0) 


0.0000 


0.292 


-0.962 


4.093 


E/Eo 


1.0000 


0.410 


0.310 


3.493 
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Table 4. Summary of Papers I-III and the present work 



paper 


disk 
geometry 


background 
feature 


magnetic 
field 


instability threshold 
and outcome ^ 


physical 
mechanism 


paper I 


razor thin 


none 


azimuthal 


Qc = 1.2 -lA 


high shear: Swing 




disks (2D) 




iPy = OO, 10, 1) 


Mcioud = a few 10'^ M© 


low shear: MJI 


paper II 


razor thin 


spiral 


azimuthal 


Mcioud = 4x lO^M© 


MJI inside 




disks (2D) 


arms 


iPy = OO, 10, 1) 




spiral arms 


paper III 


finite thickness 


none 


azimuthal 


Qc ^ 1.0 


Swing; MJI with 




disks (3D) 




{I3y = OO, 1) 


Mcioud = a few 10^ Mo 


Parker-driven boost 


this paper 


finite thickness 


none 


vertical 


Qc = 1.6 


Swing of MRI-driven 




disks (3D) 




(/?, = 100,400) 


Afcioud = a few 10^ Mq 


turbulence 
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Fig. 1. — Time evolution of maximum density 
lines) of (a) /3 = 100 and (6) /3 = 400 models, 
growth rate of the most unstable MRI mode for 



(solid lines) and density dispersion (dotted 
The dashed line in (a) represents the linear 
the parameters of non-self-gravitating model 
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Fig. 2. — Gray-scale images of the perturbed density Sp/po{0) at four time epochs from 
model A. Each frame (a)-(d) displays slices at y = —Ly/2 {left; x — z plane) and x — —Lj./2 
{right; y — z plane). The MRI "channel solution" that appears strongly at t/torh — 3 {frame 
[hj) breaks up and generates MHD turbulence at t/to^h ~ 4 {frames [c],[dj). 
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Fig. 3. — Vertical variations of horizontally averaged quantities in model A. The dimen- 
sionless, differential surface element dS is defined by dS = dxdy / {L^Ly) . (a) The averaged 
density distribution remains close to the initial profile even when turbulence is fully devel- 
oped, [h] In the saturated state, (3 is reduced to ~ 1 — 10 at \z\/H ^ 2. In (c), the magnetic 
energy density dominates the total energy density, and in (c?), the Maxwell stress dominates 
the total X — y (radial-azimuthal) stress. 
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Uo.o 



X 
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Fig. 4. — (left) Horizontal velocity and (right) magnetic field vectors at the midplane of model 
A at t/toTh = 5.6 are overlaid on the maps of midplane density {greyscale; units po(0)). 
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Fig. 5. — Amplitudes of the power spectra along the kx, ky, and k^ axes of density (left) and 
azimuthal magnetic field (right), averaged over i/^orb = 5 — 6, in model A. Along each axis, 
modes indicated by different lines are as follows: for k^, riy — riz — (solid), Uy — l,nz — 
(dotted), Uy — 0,nz — 1 (dashed), Uy — Uz — 1 (dot-dashed); for ky, Ux — riz — (solid), 
Ux — l,nz = (dotted), Ux = 0,nz = I (dashed), = = 1 (dot-dashed); for kz, 
Ux = Uy = (solid), Hx = l,ny = (dotted), Ux = 0,ny = 1 (dashed), Ux = ny = 1 
(dot-dashed), where Hx = kx/2'KLx, ny = ky/2'KLy, and = kz/2'KLz. For comparison, the 
slope of the Kolmogorov spectrum (~ k~^^/^) is indicated in each pane. 




Fig. 6. — Surface density maps projected on the x-y plane of model C. Numbers labeling 
gray-scale bars correspond to log S/Sq. A square in the right panel indicates the sector 
viewed as a three-dimensional visualization in Figure 8. 
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Fig. 7. — (a) Mean specific angular momenta inside surfaces defined by critical density pcrit 
of self-gravitating clouds in three-dimensional simulations. For each cloud, is normalized 
to Jgai (see eq. [9]), the total specific angular momentum that the cloud would retain if angular 
momentum were conserved during quasi-isotropic formation. MRI indicates model C in the 
present work, while Parker and Swing refer to magnetized model D and unmagnetized model 
C of Paper 111, respectively. Note that the magnetized clouds clearly lose specific angular 
momenta as they collapse, while the unmagnetized model does not. 
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Fig. 8. — Perspective visualization of an isodensity surface (p/po(0) = 10) and selected 
magnetic field lines in model C at t/torh = 5.5. The region fits in the projected square 
marked in Figure 6, with a vertical extent of |^| < 0.04^2. The magnetic field lines in green 
run from the right edge of the box into the clump and from the clump to the left edge, while 
blue indicates field lines that extend to the ceiling of the box. Field lines that return back 
to the clump are drawn in black. The midplane density slice is shown in colorscale at the 
bottom of the box (colorbar labels log p/ Po(0)); horizontal velocity vectors in the same plane 
are also drawn. 
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